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ABSTRACT 



Methods are developed, and Fortran 63 CODAP computer programs 
are demonstrated, to generate random numbers from the uniform, 
normal (including multivariate normal), Poisson, and exponential 
probability distributions. Various statistical tests are described 
and the results of the application of these tests to the generators 
are tabulated. A general method for generating random numbers from 
a large class of distributions is described. The methods of 
generation are optimized to provide an accurate generator while 
producing numbers at a maximum rate. The uniform generator that is 

• r 

used as a basis for the other generators is of the congruential 
type and is capable of generating 1800 numbers per second. 
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Introduction 



The need for a rapid reliable source of random numbers from 
a prescribed distribution is well recognized, and is likely to 
become even more pressing in military circles in view of the 
Department of Defense's increased interest in Operations Analysis. 
This thesis presents methods and programs, some old, some new, 
for generating random numbers from various specified distributions. 
Some statistical tests of the programs and their results are 
described. 

All methods assume a uniform random number generator is 
available. A thesis by Barron [lj has a good bibliography on 
this subject. The generator used here, and described more completely 
in Section 3, is of the mixed congruential type. While some uniform 
generators may have advantages over the one used, this one seems 
to perform very well, at the same time as being as fast as any 
demonstrated. Since this generator is used as a basis for all the 
others it should be remembered that no generator can be considered 
'perfect', especially in the continuous distribution case, since the 
computer is limited to a finite set of possible numbers. However, 
for practical purposes this inaccuracy is not important. Other 
sources of inaccuracy c an be important, however. The numbers 
generated must have the property of randomness, and must faithfully 
represent the desired distribution. These properties are measured 
by testing samples of the generated numbers. Some of the methods of 
generating numbers in this paper theoretically provide an exact 



1 



transformation from the uniform distribution to the desired 

distribution. An example of this is the half -Gaussian method 

1 

of generating normal random numbers . If we assume that the 
uniform numbers are accurate then we are led to the conclusion 
that the normal numbers are also accurate and that tests of these 
numbers would be superfluous. However, tests are performed to 
assure that the uniform numbers were so good that they did not 
bias the derived distribution. Other techniques such as Marsaglia's 
technique (see Section U) for generating normal random numbers 
are in a sense curve-fitting techniques and only provide a 
controllably good approximation to the real distribution. The 
advantage of the approximation techniques is in the far greater 
speed with which they may provide the desired numbers. The user 
who needs to draw numbers from a distribution he suspects is 
normal, with inaccurately measured mean and hypothesized variance, 
is not in need of a generator that is accurate in the sixth decimal 
place, however, he still would like to be assured that having 
assumed a distribution and its parameters he will be able to generate 
numbers with the appropriate shape and with good properties of 
randomness . 



1 The 



le phrase ' normal random numbers 1 and other like phrases 
should be read as 'numbers distributed as if coming from the normal 
distribution. ' 
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2. Testing the Generators. 

2 .1 General Discussion. 

The uniform generator chosen here has been tested extensively 
by others. The tests that have been applied include frequency 
tests, serial tests, moment tests, poker tests, gap tests, and 
many others. As mentioned in the introduction, this generator is 
as good as can be found, considering the requirement for speed. 
However, the derived distributions will be tested to overcome any 
doubts there may be about the transformation. The numbers will 
be tested mainly to measure the faithfulness with which they 
represent the derived distribution. The randomness is provided 
by the uniform generator. If the randomness is not satisfactory 
then the uniform generator must be blamed, not the transformation. 

A slower but more satisfactory method is available if the need is 
felt."*” Thus the tests used here, the moment test, the hypothesis 
tests on the mean and variance, and the Kolmogorov-Smirnov goodness 
of fit test, are not designed to detect special types of non- 
randomness, such as is detected by the poker test and other similar 
tests. 

Martin Greeriberger has written an interesting article [17] on 
this subject. He presents the results of an investigation by 
Joseph Lach [l£Q at Yale University in which Lach showed that the 
congruential method of generating uniform random numbers has a 

^See page 22* 
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predictable non-randomness in second order serial correlation. 

The lesson is clear. Having found an unders'irable feature in a 
generator, it is generally possible to modify the generator to 
eliminate the feature, however, we can be sure that we have 
introduced another aberration of some kind, even though its form 
may be hard to determine. When the user asks, "Is this generator 
good enough?", the obvious retort is "good enough for what?" 

No one generator is suitable for all applications, but the generator 
used here will be good enough for most. If the user thinks that 
this is not so in his application, he has at his resources the 
modifying methods of Marsaglia C 8 ] or Lach with the penalty of 
longer generation times. 

2.2 Moment tests 

The first four moments are calculated and are compared with 
the theoretical moments for each of the distributions. A more 
appealing statistic than the sample moment may be the unbiased 
estimator of the moment, however for large sample sizes such as 
are used here, this varies very little from the sample moment, and 
the sample moment is easier to handle in other uses- such as hypothesis 
testing. The unbiased estimator of the variance is: 




The statistic used here is the sample second moment: 




^All numerations arc on the index '4.' which runs from 1 to N, 
the sample size, unless otherwise indicated. 



a 



The first moment is the mean: 

The third and fourth moments are: » 

(£xiT J 

(2-X?) (£*;)*] 

2.3 Hypothesis tests on the mean and variance. 

The availability of large sample sizes is used in 
designing tests on the mean and variance. The cental limit 
theorem is used where possible to simplify the test procedures. 

2.3.1 Test on the mean. 

This test is applied to the normal distribution. The 
hypothesis is that the mean is zeroj, the alternate hypothesis 
is that the mean is not zero. The test is performed by 
calculating the statistic Y, where 

Y= Jn * 

The decision rule at the alpha level becomes: 

Accept the hypothesis when — Y < 

At the 10$ level this rule becomes: 

Accept the hypothesis when — l>£> t rS< lit 4*5 
The justifications for using this test in preference to the *T 1 test 
are that the variance can be assumed to be one, and that a large 
sample size is being used. 

2.3.2 Test on the variance. 

Again for the normal distribution, the hypothesis is that 
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the variance is onej the alternate hypothesis is that the 

variance is not one. The test is performed by calculating the 

statistic Z, where jz — S-(M-I) 

i|2(N-0 

s = £(x;-x ) 2 

The decision rule at the 10$ level becomes: 

Accept the hypothesis when — l.tU-S < 

2,h The Kolmogorov -Smirnov goodness of fit test. 

In Massey's discussion of this test [6] , he presents evidence 
to indicate that this test may in many circumstances be better 
than the more usual chi -squared goodness of fit test. To test 
the hypothesis that a sample of size N comes from theorized 
distribution, the cumulative step function S N (x) is formed. 

S N (x)“k/n 

where k is the number of observations less than or equal to x. 

The selection of x is arbitrary within certain limits. In this 
paper x is chosen so that there are either twenty or fifty equal 
intervals spanning the sample space. S^(x) is compared with the 
theoretical value of the cumulative distribution, F(x). The 
maximum difference d is calculated. 

d=raax|F(x)-S N (x)| 

Tables due to Smirnov [7] give certain critical points of the 
distribution for various sample sizes. For sample sizes over 
35, and at the 10$ level of significance, if 

d/n < 1.22/JTT 
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the sampled distribution is accepted as the hypothesized 
distribution. 

2.5 The chi -squared goodness of fit test. 

This test was developed by Pearson and represents the 
earliest non-parametric decision-making test in statistics. 

Partly because of the length of time the test has been in 
existence, many drawbacks to the test have been noted, however 
it still stands as a useful and much used test. For this test 
the sample space has been divided into k intervals and the number 
of sample observations in each interval is noted. 




where cj* is the observed number of sample observations in each 
interval, and m* is the expected number. Pearson showed that 
X^has the chi -squared distribution with k-1 degrees of freedom. 

The decision rule at the alpha level becomes: if X*— 

accept the hypothesis that the distribtuion is as postulated. There 

are several problems associated with the application of this test. 

How should the interval size be determined? How many intervals 
should there be? Mann and Wald [ill] studied this problem and 
formulated a criterion for the selection of k, Williams [l5j notes 
that this criterion is not particularly sensitive to even a reduction 
by a factor of two in the number of intervals. The use of equal 
probability intervals vice equal length intervals is also recommended. 
However, the basis for this recommendation is unclear. It is 
agreed that very low probability intervals such as would occur 
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in the tails in an equal interval length division of the normal 
density function should be avoided. The test is applied here 
using equal length intervals. Low probability intervals are 
avoided by 'pooling' several intervals until the probability is 
of the same order of magnitude as in the other intervals. 

2 .6 Scatter diagrams 

The best type of test to apply to the generator initially 
is some type of scatter diagram. The scatter diagram can often 
immediately give an intuitive idea as to whether the generator is 
behaving properly. In fact' the scatter diagram can be a very 
powerful tool for rejecting a generator- more sophisticated 
techniques are needed to accept the generator, however* The 
scatter diagram that was used here was constructed by plotting 
the first number generated versus the second, the third versus 
the fourth, and so on. This type of plot will also enable us to 
look for correlations similar to those found by Lach £l8j and 
discussed further in Section 2.1. 
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3. The Uniform Distribution. 

3.1 Disbribution characteristics. 

It is desired to generate numbers such that: 







x £ o 
o<-x < i 



The first four moments are : Ml = 1/2 : M2 = 1 /12 ; M3 = 0 
Mli = 1/80. The uniform distribution is usually specified 
terms of its interval - uniform on the interval from zero 
( denoted here by U (0,1)). In the general case U (A,B), 
transformation from U (0,1) is: 



in 

to one 
a simple 



URAB » (UR01) (B-A) + A 

where UR01 is the number provided by the U (0,1) generator, 
and URAB is the random number uniform on the interval (A,B). 

3.2 Methods of generation. 

ft 

3.2.1 Many techniques have been used over the years. For some 
particular applications such a method as table look-up may be 
suitable. However for our purposes what is desired is a rapid, 
•accurate' method for the computer to produce a practically 
inexhaustible supply of numbers. 

3.2.2 An early computational scheme was called the mid-square 
method. In this procedure two starting values, say A1 and A2, 
are multiplied together? the middle set of bits (usually 2h) are 
extracted as the third random number A3j then A2 and A3 are 
multiplied together and the algorithm is continued in a similar 
fashion. This method has performed well in many tests but 
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unfortunately degenerates to all zeros in a relatively short 
time. 

3.2.3 An improved computational scheme was tested extensively 
by Hull and Dobell [3] . This method forms a series of numbers 
Ai, where 

A^ + -]_ = BA^ + C modulus 2 

B and C are constants to be selected. 

3.2 .U A special form of this generator where C = 0 is the 
generally recommended form. C is set equal to zero because it 
does not improve the characteristics of the generator and it adds 
to computation time. The selection of B and the starting number X 
has been the subject of considerable study. Barron and others 
have shown that: 

= (2^ + 3 )Xr modulus 2 ^ 

where X is either 1, or 2^® -1, or any number naturally generated in 

the sequence, is an excellent generator. Since this generator had 

been tested extensively previously no attempt was made to test it 

rigorsly although several interesting characteristics were noted. 

Some of these are included here. The generator was run down through 

7 8 

the first 10 to 10° numbers. A number at the end of this sequence 
was extracted for possible alternate use as a starting number. It 
can be found in Appendix I. A graph of the mean of the first 10000 
times i numbers, where i runs from 1 to 100 is plotted against the 
index i. This plot is compared with curves of = K^q^/ lOOOOi + 0.50000 
versus i. The more our data plots between these curves the better. 
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We would expect it to be between the curves 90 % of the time. 

As the following graph shows our generator does not quite 
live up to this expectation. A scatter diagram of the type 
suggested by Lach [l8] was plotted. The graph on page 13 
consists of 3000 points constructed as described in Section 2.6. 

3.2.5 George Marsaglia and M. Donald MacLaren [8] at Boeing 
Scientific Research Laboratories suggest that the combination of 
two generators will produce a superior random number of generator. 
They have tested: 

Ar + ]_ = (2^ + 3)Ar modulus 2 
and Br+]_ = (2^ + l)B^ modulus 2-^ 

In essence they have used one generator to select numbers from the 
other. This generator seems to provide an improvement in some 
local randomness properties. However, the penalty for the improve- 
ment is doubling the time of generation. Marsaglia and MacLaren 
also noted that the method of table look-up may once more become 
feasible. In the case of the CDC 160U this method is not practical.. 
However, in a parallel program c ompu ter a method using a pair of 
generators may well be advantageous. The generators continually 
fill up the bottom of a short table in memory as the main program 
uses numbers from the top. The size of the table is chosen to 
ensure that the program never uses all the numbers in the table 
and thus the effective generation time will be just the load cycle 
time. 

3.2.6 The CODAP, Fortran 63 , machine language program for the 
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Cumulative Mean Of the First 1* Samples 
Of The Uniform Random Number Generator 
(with 1000 numbers per sample) 
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Scatter Diagram 

For The Uniform (0,1) Random Number Generator 
(Scaler 0.2 units per inch) 
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uniform generator used as a basis for this thesis is in Appendix 
I. The expected time of generation per number, as calculated over 
several samples of varying size was found to be 552 microseconds 
per number. This amounts to producing l8ll numbers per second. 
However, the generator is theoretically much faster than that. The 
time per number as calculated from times in the Control Data 

z •- 

Corporation specifications for the computer is 121 microseconds. 
Measured times, depending on the context and the timing mechanism 
varied from a minimum of 370 to a maximum of 700 microseconds. 



u. 



The Normal Disbribution (Univariate Case). 



U.l Distribution characteristics . 



It is desired to generate numbers such that the 



density function will be 




**r(-k (*£*)') 



where u is the mean of the distribution and a 2 - is the variance. For 
the basic case the mean is taken to be zero and the variance to be 
one. The first four moments are Ml = 0, M2 = 1, M3 = 0, Ml; = 3* 

If the desired distribution is to have a mean other than zero, 
say u, and a variance other than one, say V, then the following 
transformation applied to the numbers generated by the N(0,l) 
generator developed here, represented by RN01, will, produce a 
number, RNUV, With the desired characteristics. 



In this paper the normal distribution is treated in three 
separate sections. The univariate case is developed first, then 
the bivariate, and finally a general multivariate case is demonstrated. 
The main purpose of this separate treatment is to allow a more 
efficient handling of the more commonly used univariate and 
bivariate cases. A general n-dimensional normal random number 
generator would be much slower, when used for n equals one, than 
the univariate generator demonstrated in Section U. The test 
procedures for each generator are also different. 

U.2 Methods of generation. 

U.2.1 The normal distribution is one of the most used and 



RNUV = (RN01) (V) + U 
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tabulated distributions. As with the uniform distribution several 
procedures have been used to produce random numbers from it* The 
methods discussed here are those that are most adaptable to use 
on a computer. 

U.2,2 The most common procedure has been to take the sum of K 
uniform random numbers. The central limit theorem shows that this 
(with the mean subtracted, and divided by the standard deviation) 
approaches the normal as K gets large, Vaa tested a generator 
using the sum of twelve uniform random numbers. This approximation 
has the disadvantage of being truncated at plus and minus six. 

Even more important a factor is the time required to generate these 
numbers . It is hoped that a more exact and faster method can be 
found, 

U.2.2 The so-called half -Gaussian method [l]] provides a theoretically 
exact transformation from the uniform to the normal. However in an 
attempt to reduce time some fairly drastic approximations have been 
made. These approximations should not affect randomness, however, 
and should only be a factor in accuracy beyond the fifth decimal 
place. The following flow chart is the basis for the routine. R(J) 
is a uniform random number. HGRN is the normal number. 
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The routine first generates the positive half of a normal distribution 
then adds a random sign selected by another uniform number. This 
program makes no external calls to the log function but rather 
uses the following series approximations 

T = (R-1)/(R+1) 

LN(R)= 2(T+t3/3+ T 7 /7+ T 9 / 9 ) 

The CODAP function sub-program is Appendix II. 

b.2.3 An excellent approximation technique has been developed by 
Marsaglia and Bray £ 9 ] . Marsaglia has developed several similar 
techniques but the one proposed here seems optimum in terms of time 
required per number and the storage space required. The method 
involves selecting one of four functions of varying complexities to 
produce the random number. One function is very simple and fast and 
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is used 86$ of the time; the next function is also fast and is 
used 11$ of the time. The remaining three percent of the time 
much more complex functions are used, however, due to the rapidity 
with which 97$ of the numbers are formed, the overall expected 
generation time per number is relatively low. The program outline 
is as follows (MSRN is the desired random number): 

1. 86.38$ of the time, set 

MSRN = 2(R(J) + R(J+1) + R(J+2) - 1.5] 

2. 11.07$ of the time, set 

MSRN = 1.5 (r(J) + R(J+1) - ^ 

3. 2.28002039$ of the time form pairs (X,Y) such that 

X = 6R(J) - 3 and 

Y = 0.358 R(J+1) 

until Y*G3(X); then set MSRN = X. G3(X) is defined by: 
17.U9731196exp(-xV2)-U.73570326(3-x i )-2.l57875UU(1.5-lx|) for IxKl 
17.U9731196exp(-xV2)-2.36785l63(3-lxl) Z -2.l57875UU(1.5-lx| ) for l<x<1.5 

17 .U973H96exp(-x 2 /2)-2. 36785163(3-1x0 for 1.5<x<3 

U. 0.26997961$ of the time form pairs (X,Y) until either X 
or Y is greater than three, then let that one equal MSRN. For 
RM(J) uniform on the interval (-1,1) and such that if RM(J) X + RM(j+l)*Sl, 
then let Z = RM(J)*+ RM(J+lf 
and X = RM(J)[ f 9-2LN(Z)j/(Z)] 

Y = RM(J+1)[ (9-2LN(Z)j/(Z)] 

The CODAP function sub-program is Appendix III. 
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ii.3 Selection of generators. 

A complete table of results of the tests on the normal 
generators is in Appendix IV. Partial results are presented below. 

HALF-GASSIAN TECHNIQUE MARSAGLIA TECHNIQUE 

Average observed 362$ microseconds 1503 microseconds 

generation time 
per number 



Sample range 


1-1000 


1-10000 


10001-20000 


1-1000 


1-10000 


10001- 

20000 


Sample size 


1000 


10000 


10000 


10000 


10000 


10000 


Mean 

(theor. = 0) 


0.0U 


-0.01 


0.00 


0.00 


-0.02 


-0.01 


Variance 
(theor. 14 l) 


1.06 


l.Oli 


1.01 


.99 


1.00 


.99 


3rd Mom. 
(theor. = 0) 


- 0.05 


-0.02 


0.01 


.01 


-0.0U 


-0.07 


Uth Mom. 
(theor. =3) 


3.17 


3.20 


3.12 


2.77 


2.97 


2.93 


JN X 

(K <o5= 1.6U) 


1.32 


- 0.96 


0.30 


0.13 


-2.U7 


-1.06 


S-(N-l) 

i(N-l) 


1.30 


2.57 


0.95 


-0.1U 


-0.11 


-0.80 


^max/'* 


0.0270 


0.0073 


0.0033 


0.0082 


0.0095 


0 . 00 U 9 


1.22 /JW 


0.0386 


0.0122 


0.0122 


0.0386 


0.0122 


0.0122 



The samples generated by the half -Gaussian technique passed the hypothesis 
on the mean, at the 10$ level, 80$ of the time; on the variance at 
the 10 $ level, 30$ of the time. At the 5$ level the test on the 
mean was passed 90$ of the time, the test on the variance was passed 
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5 0 $ of the time. The Marsaglia technique generated samples 
that, at the 10$ level, passed the test on the mean 80$ of 
the time, and on the variance 90$ of the time. The Marsaglia 
technique consistently passed the Kolmogorov-Smirnov test but 
appeared a little heavy in the ’tails' none the less. For a 
further examination of this see Addendum 1. The half-Gaussian 
technique failed the Kolmogorov-Smirnov test once but appeared 
better behaved in the tails. The graph on page 21 is a scatter 
diagram for the numbers produced by the Marsaglia technique. 

The graph consists of U£00 points. Marsaglia 's technique produced 
better results for the first four moments. The decisive factor 
that leads to the selection of one of the generators is the time 
to generate each number. The Marsaglia technique is more than 
twice as fast as the half-Gaussian technique, is the one used in 
further developments, and is the one recommended for general use. 
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Scatter Diagram 

For The Normal Random Number Generator 
(Scale; 1.0 units per inch) 
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5. The Normal Disbribution (Bivariate Case) 

5.1 Disbribution characteristics. 

It is desired to generate pairs of numbers >X Z ;, ) 

such that the two-dimensional random variable (Xq_, Xg) has the 
joint density function 



for all (x^^). The constants are u^ and i^, the means ; 0 \ (> 0), CT X 
(>0 ) , the standard deviations; and (> (-l^il) , the correlation 
coefficient. Thus the requirement for a random vector must be 
accompanied by the specification of the mean Vector (u^i^), the 
variances (the squares of the standard deviation), and the 
correlation coefficient. Another equivalent form of the input 
would be the mean vector and the covariance matrix. This last 
form of input will be used in the general multivariate case. 

The distribution is specified in matrix notation as follows: 



f(x) . - ^ecp(-i(Y-U)'R(r-u)) 

(2n ) *■ 

where f(X) is the joint density function of the x's; p is the 
dimension-in this case two; U is the mean vector; and R is the 
inverse of the convariance matrix. Thus iRl'^is the square root 
of the determinant of the inverse of the covariance matrix; (Y-U) 
'R(Y-U) is a quadratic form, where (Y-U)' is a one by p vector, 

R is a p by p matrix, and (Y-U) is a p by one vector. 
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5.2 Method of generation. 



If R(J) is a normal random number, the random bivariate 



vector (VN1, VN2) is formed as follows: 
VN1 = ( cr x )R(J )+ U| 



VN2 = ( 07 )R(J ) + ( 0^)R(J+1)( 



The source of the normal random numbers is the Marsaglia routine 
described in the previous section. The generator is Appendix V. 

5.3 Testing the generator. 

First the maximum likelihood estimators of the mean vector 
and the covariance matrix are formed [12] . The maximum likelihood 
estimators for the parameters are: 



The distribution of the mean when the covariance matrix is unknown 
was shown by Hotelling to be a multivariate analogue of the t-test 
and is called the generalized T statistic. However, the covariance 
matrix is known and once again we can use a more powerful test. 




<© a - £ 



;-*) J 

($j 7 - xfSf 1 (6,2 . 



H o : u = (U 10 U 20 )' H A : u* (U 10 U 2q ) 



Construct the statistic H such that: 



H = N(x!-u 10 , x 2 -u 20 )C" 1 (x 1 -u 10 , x 2 -U2 0 )« 
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where C is the given covariance matrix. If H < X * 0*)we accept the 
null hypothesis. The test for a hypothesized mean vector 
U = (0 , 0) ' , and a covariance matrix with variances one and 
correlation coefficient £ , reduces to calculating H such that: 



H = 






X, 



- L ( X?-- + *2* ) , for cr, rcr t -l , 

\~e z 1 

for sample sizes of 1000 the test will be made for^ = 0 . 25 , 0.5, 

0.75. At the 10$ level if H — U.6l accept the hypothesis. In 
order to test whether the covariance matrix is a given matrix, 
the test as outlined in Anderson C12J could be used, but before 
this a digression into the philosophy of testing is in order. 

The purpose of these hypothesis tests is in general to render a 
judgement as to whether a sample of data from some type of 
experiment is distributed according to a certain distribution 
with certain parameters. This type of test was applicable to our 
uniform random number generator, but its use seems extraneous for 
the distributions derived from this generator. The trans- 
fer .lions from the uniform to the normal to bivariate normal are 
either theoretically exact or controllably close to being exact. 

What is really needed is a method to ensure that any 'inaccuracies 

there may have been in the basic generator are not somehow 

amplified in the transformation so as to bias the derived distribution. 



To this end sophisticated testing procedures do not seem to be 
in order. Instead the testing will be restricted to calculation 
of estimators, and some goodness of fit tests. The testing plan 
for the bivariate generator is to test 10 sets of 1000 each for 
each of several correlation coefficients. The generator was also 
tested for some odd combinations of parameters- such as U= (100.0 
-100.0), ofsO. 5, erf- 25.0, (0= -0.8. 

With 0.7$ as the correlation coefficient, the fourth set showed a 
disagreeably large difference in the maximum likelihood estimators. 
The H statistic was close to the critical value at the 10$ level. 
The sample passed the test in 90$ of the cases. With 0.50 as the 
correlation coefficient the sample passed the H test 90$ of the 
time. Sample set four once again had rather low covariance 
estimators and once again had 3»?6 as the value for H (compared 
with h.6l for the critical value). Sample set six was again the 
culprit in not passing the H test. A similar result was noted 
in the set using 0.25 as a correlation coefficient. This suggests 
the advisability of further testing of the Marsaglia generator. in 
these ranges. Complete test results are in Appendix VI. The 
generator performs as expected and is recommended for use. 
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6. The Normal Disbribution (Multivariate Case). 

6.1 Distribution characteristics. 

As is noted in Section 6.1 the desired distribution is 
f(X) = f(x,, Xi ,...,x N ) = — f / x ecp(- 
where R is the inverse of the covariance matrix. 

6.2 Method of generation. 

As shown by Wold [l6] , the method is based on a triangularization 
of the covariance matrix such that if: 



Cll C12 ... C1N 
C21 C22 • • • C2N 

* t 




‘Pll 0- • 0 • 
P21 P22 •* 0 

• * v ' 




"Pll P21 *■- PN1 
0 P22 '* • PN2 


* • 




V % » X 




• 

> • * 


* ... 
.CN1 CN2 ' ' CNrt 




% . ♦ . 
PN1 PN2 PNN 


< 


* v • 

0 - 0 PNH 



then the N-dimensional random vector Z may be formed as follows. 

The RN(l) are the normal random numbers generated by the Marsaglia 
technique. 

Z(l) = Pllx RN (l) 

Z (2 ) = P21 x RN(1) + P22 x RN(2) 

Z (3 ) = P31 x RN (l) + P32 x RN(2) + P33 x RN(3) 

Z(N) = PN1 x RN(l) + PN2 x RN(2) + ... + PNN x RN(N). 

The method as programmed, assumes all the means are zero, but simple 

addition of the mean when required will remedy this. The matrix 
triangularization is based on a symmetric, positive definite or 
positive semi -definite matrix C. Thus any pair of the random 
variables can have a partial correlation coefficient of one. The 
routine will set all elements of the vector that are dependent equal 
to zero. This procedure was selected since in order to relate the 
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variables properly requires an inordinately extra amount of 
computer time- time that even when not needed adds to execution 
time. If the user desires some variable to be a linear trans- 
formation of some others, then he only needs to keep track of 
which variables these are, and where the program sets the vector 
element equal to zero, substitute the appropriate linear 
combination of the corresponding independent elements of the 
vector. The routine also checks and where the dependence is very 
close to one will assume a correlation of one, and proceed as 
noted above. This procedure is required to prevent division by 
numbers very close to zero. The following example will clarify 
the above explanation. Suppose it is desired to generate random 
vectors from a distribution with mean vector one and covariance 
matrix C, where: 



2 1 3 ' 



C 



l U 5 
3 5 8 



Thus column three is a linear combination, in fact the sum^of 
columns one and two. This means that the third variable is not 
an independent variable. The triangulation routine will produce a 
matrix P of the form: 
a o o' 
b d o 
c e o 
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As expected column three is identically zero. Thus if the first 
two normal (0,1) numbers generated are denoted by R1 and R2, the 
generated vector will be of the form: 

Z = (aRl, bRl + dR2, 0) . 

Since it has been determined that the third variable is the sum 
of the first two and also that we desire all the means to be one 
then the desired vector is of the form: 

Z = (aRl + 1, bRl + dR2 + 1, aRl + bRl + dR2 + l) 

The generator is Appendix VII. 

6.3 Testing the generator. 

The maximum likelihood estimators of the mean vector and 
the covariance matrix are formed as follows: 

BM1(J) = Tvf S Y( r < J J for 

C(I,J) = ^[1 (Y(r,K)*r(T,X)l] for r,T=>, 

N is the dimension of the covariance matrix, M is the number of 
sample vectors, the Y(I,J) are the vector elements, the BMl(j) 
are the elements of the mean vector estimates, and the C(I,J) are 
the elements of the covariance matrix estimate. The results of 
the tests for various covariance matrices are listed in Appendix VIII. 
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7. The Poisson Distribution. 

7.1 Distribution characteristics. 

It is desired to generate numbers such that: 

P fx “ a] = e “ m m a for m>o;a = al, ... 

~af~ 

Thus the Poisson distribution is a discrete distribution for 

integer values of a, and is characterized by the parameter m. 

2 

The first four moments are Ml ** M2 «* M3 = m, and MU = 3^ + m* 

7.2 Method of generation. 

This method is due to Kahn cm and is a theoretically 
exact transformation. The flow chart of the routine is shown below: 




POIS is the generated Poisson random number, m is the parameter r E 
is the irrational number 2*71828.,^ 

7*3 Testing the generator* 

The first four moments are calculated for samples of $000 
or 1000 with the parameter m taking on various values. The 
Kolmogorov-Smirnov goodness of fit test is applied to some of the 
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samples. As Tate and Clelland [l9j have stated the test is 
applicable to discrete distributions with neglible changes in 
significance for large sample sizes. The generator performed 
well and consistently passed the test. The test results are 
in Appendix X, while the generator itself is in Appendix IX. 
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8 . The Exponential Disbribution. 

8.1 Distribution characteristics . 

It is desired to generate random numbers such that the 
density function is: f(a) = e~ a . 

The first four moments are Ml = M2 = 1, M3 = 2, Ml: ® 9 * The 
distribution is often specified as: 

f( a ) = X e“ Aa for ^ 

where A is the parameter of the distribution. The generator here 
takes the case where A = 1 * However the exponential distribution 
has the characteristic that if the numbers generated here, (Expl) 
for which A = 1 , are simply multiplied by the parameter desired 
for the distribution (tAMBDA), then the desired numbers are generated. 
EXPL = EX PI * LAMBDA 

8.2 Method of generation. 

This method is a theoretically exact procedure due to 
Marsaglia fl3] • A more obvious method would be to take the 
integral transformation- the negative logarithm of a uniform 
random number. However this method is slower on most computers 

r 

than the one demonstrated here. Let C = l/(e-l), and let the 
random variable N take on the values l,2,3,ii,... with probabilities 
c, c/2! , 0/3 I , . . . Then let the random variable M take values 0,1, 

2 *5 

2,3,... with probabilities 1/ce, 1/ce , 1/ce- 5 ,... Then we form the 
desired random number: 

EXEL = M+MIN (Ul, U2,...,UN) 

The CODAP generating routine is Appendix XI. The sample moments and 
the results of the goodness of fit tests are in Appendix XII. 
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9. A General Disbribution. 

The method of obtaining random numbers described in this 
section is applicable only to a restricted set of distributions. 
However, the method is applicable to a type of distribution that 
frequently occurs in model building and war gaming.. The method 
is examined by use of an example. A discussion of how and when 
to use this method for other distributions is included. 

It is desired to produce numbers from the density function 
f(x) such that: 

f(x) = 2-2x for 0 - x-^lI 

The method is based on drawing uniform numbers in pairs, normalizing 
the scale of the numbers, and testing to see if the point formed 
by the pair lies under the curve described by the density function. 
If it does, the x coordinate of the point is taken as the random 
number; if it does not, another pair of uniform numbers are drawn 
and the process is continued. The flow chart for f(x) = 2-2x is: 




The U(J) are the uniform (0,1) random numbers and RN is the random 
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number from the distribution f(X). Thus the method could be 
applied with theoretical exactness to any bounded continuous 
distribution. Another commonly used density is g(x) where 
g(x) = nx 11 for 0<x<l, n - 1 

This density function is bounded and continuous and a member of 
the class to which this method is applicable. In many applications 
the user must make a judgement about the amount of use a generator 
is going to get. If it is intended for heavy use it may be well 
to explore the literature for, or to design, a more efficient 
generator than the type described in this section. However, if the 
generator is to only be used a limited amount, or if only a 
restricted amount of resources are available, this generator is 
easily programmed and will be eminently satisfactory in a wide 
variety of cases. The efficiency of the generator may be defined 
as the reciprocal of the expected number of iterations needed to 
produce a point under the curve. 





Figure la clearly has an efficiency of 1/2, figure lb has an 
efficiency of less than 1/2, and as n gets large the efficiency 
decreases rapidly. 
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A distribution such as figure lc where h(x) goes to zero at some 
large value of x may vary well have such a low efficiency as to 
make this generator unsuitable. If h(x) only approaches zero as 
x becomes large, and therefore x is not bounded, then another 
generator should be used. However, if one cannot be found and 
the user is not too concerned with the tail of the distribution, 
we can easily adapt this generator. Suppose we wished to generate 
numbers for the Weibull distribution, which is a three parameter 
distribution used widely in reliability theory, then the following 
procedure might be followed: 



Since the distribution may take any of the forms in figure 2, we 
must first limit the development to those parameter combinations 
that are bounded at x = 0. Since the user is also often concerned 
with the behaviour in part of the tail it must first be decided if 
the distribution could simply be truncated at some point. Even if 
this is acceptable an efficiency of .00001 can easily be envisioned. 
If it is unacceptable, the tail beyond some point could be closely 
approximated by the exponential distribution. 




2.0 



t 
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10. Summary and Conclusions. 

The analyst who desires to use any of the generators 
demonstrated here is once again warned that although the generators 
are recommended for general use they all have aberrations of some 
kind. The general form of these aberrations is noted where 
known. If a user suspects from analysis of his results that the 
aberration in the generator is influencing his results then it is 
recommended that he first modify the uniform generator by one of 
the methods suggested in section two. The chance of this occuring 
is remote and other parts of the model should be checked carefully. 

The programs demonstrated here are very fast. The problem 
of speed is stressed here and is a major decision criterion in the 
selection of generators. In some applications the speed may not 
be as important a factor, and in this case the type of generator 
discussed in section nine may be very cost effective in that 
little investment in programming and testing is required. The 
problem of measuring the speed of generation of a number is not as 
simple as it may appear. The generation time depends on the program 
using the generator and also on the method used to time the routine. 
The time to generate a number based on the Control Data Corporation 
specifications for the 16QU computer theoretically should be 121 
microseconds. The observed generation times vary from 370 to 700 
microseconds. 

The tests used here are statistically sound, but the meaning 

of the results of several tests of the same generated sample could 
bear further investigation. 
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Some very interesting new methods for generating random 
numbers from various distributions have been developed by 
H. Rubin. Rubin's work is soon to appear as a Stanford Applied 
Statistics Laboratory Technical Report* It would be of interest 
to compare his methods with the technique used above. 
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Appendix I. 

The Fortran 63 CODAP function subprogram for generating uniform 
(0,1) random numbers. 

The program is called by using a variable, 'UNIFORM(DUMMY) 1 • The 
argument 'DUMMY' is not used. For example: A = UNIFORM ( DUMMY )+ B. 

No common or dimension entries are required. 

The observed average length of the program is 552 microseconds. 

The starting number may be change to 'R OCT 5U02033U50?27lj22 ' 
if it is desired to enter the generator near the middle of the 



period. 





IDENT 


UNIFORM 




ENTRY 


UNIFORM 


Ml 


OCT 


Uooooooooooooooo 


M2 


OCT 


2000000000000000 


R 


OCT 


1777777777777777 


UNIFORM 


SLJ 






LDA 
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INA 


1 




SAU 


EXIT 




ENQ 
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LDA 


R 
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+19 




SCL 


Ml 




ADD 


R 




ADD 


R 




ADD 


R 




STA 


R 




ARS 


11 




ADD 


M2 




FAD 


M2 


EXIT 


SLJ 


•K-* 




END 





39 



Appendix II. 



The half -Gaussian technique for producing normal random numbers. 
The number is generated by the use of the expression 'GSRN(BLK)'. 
The number is returned as GSRN. The argument .'BUC' is not used. 
No external calls are made. 

The observed average length of time to generate one number is 



3625 microseconds. 



QSRN 



.190 
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GSRN 


ENTRY 


GSRN 
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LDA 
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RTJ 
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STA 


B1 


FAD 
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STA 


TEMP 


LDA 
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FDV 


TEMP 


STA 


X 


FMU 


X 
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FMU 
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X2 


FAD 
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X2 
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X2 
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=02 001 U 00000000000 
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STA 
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UNIFORM 


STA 


B1 
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=02 001 U 00000000000 


STA 
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STA 
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STA 
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THS 
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THS 
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STA' 
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X2 
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X3 


BSS 


1 


x5 


BSS 


1 


X7 


BSS 


1 


X9 


BSS 


1 




END 
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Appendix III. 

The number is generated by use of the variable 'MARS(ZQ)'. The 
argument 1 ZQ * is not used. For example: 'A = MARS(ZQ)/9. 

External calls are made to LOGF, SQRTF, and EXPF. 

No common dimension entries are required. 

The observed average length of time to generate one number is 
1503 microseconds. 



UNIFORM 



R 

ZMARS 



G1 



IDENT 
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ENTRY 


ZMARS- 
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STA 
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STA 
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FSB 
STA 
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SCM 
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STA 
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CALL 
FMU 
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LDA 
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SLJ 
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FSB 
FMU 
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STA 
LDA. 
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FMU 

FAD 

SLJ 
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DEC 
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DEC 
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END 
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Appendix IV. 

Test results for Marsaglia normal generator. 

The test consists of ten consecutive samples with 10000 numbers in each. 



MEAN 


M2 


M3 


Ml* 


JNX 


S-(N-l) 
'I z.(n-0 


d/n 


.0 


1.0 


.0 


3.0 


-1.61*5 
to 1.61*5 


-1.61*5 
to 1.61*5 




-.021*71* 


.9981*1* 


-.01*1*35 


2.96870 


-2 .1*737 


-.1101* 


0.0095 


-.01062 


.98869 


-.071*87 


2.92989 


-1.0625 


-.7996 


0.001*9 


-.00561* 


.99388 


-.07973 


2.93538 


- .561*0 


-.1*330 


0.0027 


-.00873 


1.001*52 


-.09225 


2.99577 


- .8729 


.3196 


0.0060 


-.00921 


.9991*5 


-.01*771* 


3.03286 


- .9209 


-.0392 


0.001*0 


-.00683 


.99709 


-.10706 


3.081*87 


- .6831* 


-.2057 


0.001*8 


.00028 


.99963 


-.01*737 


3.03271* 


.0277 


-.0261* 


0.0030 


.00732 


1.01508 


-.01*782 


3.161*72 


.7317 


1.0661* 


0.001*1 


-.001*10 


1.021*51* 


-.05629 


3.09820 


- .1*101* 


1.7352 


0.0050 


-.01732 


.97325 


-.05067 


2.79351 


-1.7319 


-1.8911 


0.0050 
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Appendix IV. (Cont'd) 

Test results for half -Gaussian method. 

The test consists of ten consecutive samples with 10000 numbers in each. 



MEAN 
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M3 
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S-(M-l3 
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-1.61*5 to 
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OO 
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<5 

O 

• 

1 
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-.1*71*7 
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3.29951 
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-.00372 


1.02968 
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Appendix V. 

The bivariate normal random number generator (using Marsaglia's 
technique) . 

The number is generated by a call 'CALL AKHN(VN1, VN2)'. The 
bivariate vector is returned as the arguments VN1, VN2. 

Subroutine AKHN has the following arguments in common: SIGISQ, 
SIG2SQ, RHO, Ul, U2. The SIGISQ are the desired variances, RHO 
is the correlation coefficient, and Ul and U2 are the desired 
means. Thus besides reading in these values in the main program 
they must be communicated to the subroutine by a common statement 
such as: 'COMMON/SIGlSQ/SIGlSQ/SIG2SQ/SIG2SQ/RHO/RHO/tJl/Ul/U2/U2 ' . 
External calls are made to IX)GF, SQRTF, and EXPF. 

The observed average generation time for a pair of numbers is 
6 H 80 microseconds. 
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Appendix VI. 

Test results for bivariate normal generator. 

The test consists of ten consecutive samples of 1000 numbers each for 
each of three different correlation coefficients. 
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Appendix VII 

The FORTRAN programs to produce multivariate normal random vectors. 

The first entry must be 'CALL TRIANG!. This produces the 
triangularized matrix, and allows repetitive calls to 'CALL MULTN(Z)'. 
The argument ' Z ' is the starting address of the random vector. 

Several other entries are necessary. The dimension of the 
covariance matrix (and the desired vectors) is set equal to 'NR'. 

The desired covariance matrix is stored as a matrix called 'C ' . 

A common statement with 'NR' and 'C', with 'C' appropriately 
dimensioned, is included. ' Z ' is also dimensioned. The 
following sample program will read in a matrix 'C ' , which is a 
5 by 5 matrix punched column by column on cards. The random 
vectors are stored in a 5 by 100 array. Subroutine 'MULTN' 
and subroutine 'TRIANG' follow. Function 'ZMARS', the normal 
random number generator, from Appendix III must be added. 

Subroutine TRIANG makes external calls to SQRTF. 

PROGRAM EXAMPLE 

common/Nr/nr/c /c(5,5) /p/p (5,5) 

DIMENSION Z(5), ARRAY(100,5) 

nr=5 

READ 101, ((C(I,J),I-1,NR),J-1,NR) 

101 FORMAT (5F8.U) 

CALL TRIANG 
DO 111 J =1,100 
CALL MULTN (Z) 

DO 111 K-1,5 
111 ARRAY (J,K)=Z(K) 

STOP 

END 

For a three by three matrix subroutine TRIANG takes 10600 
microseconds, and each vector is produced in 9520 microseconds. 
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It is noted that since the programs are written in FORTRAN, 
they are not in any sense optimal. 



SUBROUTINE TRIANG 

C OMMON /NR/NR/C /C ( 3 ,3 ) /P/P( 3 ,3 ) 

NC-1 

DQ 5 I - 1, NR 
DQ 5 J - 1, NR 
5 P(I,J) = 0.0 

P(l,l)=SQRTF(C (1,1)) 

IF(P(l,l) .LE. .0001)771,9 
9 DO 15 K*2 ,NR 
15 P(K,l)=C(K,l)/P(l,l) 

18 DO 83 JH=2,NR 
Q=0.0 
NC-NC+1 
NCM*NC-1 
NCP=NC+1 
DO 50 J=1,NCM 

50 Q=Q+P(NC,J)*P(NC,J) 

x=c(jh,jh)-q 

IF (X.LE.. 000005) 777,505 
505 P(JH,JH)=SQRTF(X) 

51 DO 55 K=NCP,NR 

S=0.0 

DO 53 JQ=1,NCM 
53 S=S+P(K,JQ)*P(NC,JQ) 

55 P(K,NC)-(C(K,NC)-S)/P(NC,NC) 

83 CONTINUE 
GO TO 666 
771 DO 773 L«1,NR 
773 P(L,l)=0.0 
GO TO 18 

777 DO 780 L=JH,NR 
780 P(L,JH)-0.0 
GO TO 83 
666 CONTINUE 
RETURN 
END 
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SUBROUTINE MULTN(Z) 

C OMMON/NR AR/P/P( 3,3) 
DIMENSION Z(10),RZM(10) 
DO 220 J-l, NR 
Z(J)«0.0 

220 RZM(J)-ZMARS(DUM) 

DO 230 L-1,NR 
DO 230 M-1,NR 

230 Z(L)-Z(L)+P(L,M)*RZM(M) 
RETURN 
END 
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Appendix VIII 

The tests are for a sample size of 100 vectors. 

The first matrix tested was a 10 by 10 identity matrix. The 



maximum likelihood estimate of the covariance matrix was: 
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Appendix IX. 

The Poisson distributed random number generator. 

The numbers are generated by an initial use of the variable 
'NPOISSET (mean) ' , which initializes the generator for the desired 
value of the njean. This is followed by calls to 'NPOIS(DUM)' to 
actually produce the random numbers. The argument 'DUM' is 
not used. For example for a desired mean of 3.0: 

X(l) = NP0ISSET(3.0) 



X(J) * NPOIS(DUM) 
External calls are made to .EXPF, 



The observed average generation time per number was found to be 
approximately representable by: TIME = 600 + 630(MEAN). 
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Appendix X. 

Test results for the Poisson generator. 



The test consists of ten consecutive samples of either 1000 or 
$000 numbers each-for various parameter values. 
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1.92 2 0 


1.8618 


1.5925 


10.5851* 


.02 90 








2.0120 


2.2381 


2.6755 


18.2820 


.011*7 








1.9990 


2.0110 


2.0572 


13.8852 


.0073 






theor. 


2.0 


2.0 


2.0 


11 *. 0 


.0386(10$) 
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PARAMETER 


Ml 


M2 


M3 


MU 


d/n Sample 












Size 


U.o 


3.9578 


3.8668 


3.533U 


U7.2U82 


.013 5000 




U.021U 


U.1002 


U.1851 


5U.9297 


.008 


• 


U.020U 


U.068U 


3.6521 


51.6736 


.012 




U.OU22 


U.0116 


U.082U 


51.6602 


.010 




3.9770 


U.13U9 


U.63U2 


59.U956 


.010 




3.9880 


U.0679 


U.U8U7 


58.0892 


.006 




3.9920 


U.066U 


U.7U96 


58.U535 


.005 




U.0130 


3.8916 


3.1U39 


U6.6069 


.01U 




U.0318 


U.0536 


U.17U0 


53.8986 


.008 




3.9806 


3.9U26 


3.91UU 


50.7897 


.005 


Theor. 


U.o 


U.O 


U.o 


52.0 


.017 (1056) 


5.0 


U.93UO 


U.97U2 


5.3U30 


78.80U5 


.016 1000 




5.0U62 


5.0975 


U.9371 


80.2735 


.01U 




5.0108 


U.93U5 


U.6862 


7U.9158 


.005 




U.9912 


5.0533 


5.1257 


80.938U 


.00U 




5.0038 


5.0U12 


5.7602 


87.0035 


.003 




5.0126 


5.0367 


5.339U 


83.8973 


.005 




5.03U6 


U.98UU 


U.8206 


80.172U 


.010 




5.0016 


5.01U2 


5.1778 


77.76U8 


.006 




5.02 1U 


U.9U19 


U.7638 


75.7139 


.006 




U.9598 


5.01U0 


U.8610 


77.982U 


.009 


Theor. 


5.o 


5.0 


5.0 


80.0 


.039 {10%) 


10.0 


9.9 760 


9.6751 


10.0922 


276.3385 


.013 1000 




9.8280 


9.5560 


9.8887 


280.6076 


.030 




9.9070 


10.U328 


12 . 0257 


322.7622 


.027 




10.15U0 


9.8UU1 


9.U111 


287.9662 


.019 




10.1090 


10.8179 


10.537U 


3U7.2039 


.021 




9.9710 


9.9922 


8.U5U8 


305.978U 


.012 




10.0520 


10.2676 


11.U099 


298.631U 


.029 




10.1500 


10.3258 


11.2U6U 


307.3533 


.020 




9.99UO 


9.U59U 


7.103U 


265.9739 


.017 




10.03U0 


9.6UU5 


8.7856 


283.3U59 


.011 


Theor. 


10.0 


10.0 


10.0 


310.0 


.39 (1056) 



Generation 

Time 

(microsecs) 

3128 



3735 



6935 
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PARAMETER 


MI 


M2 


M3 


Ml* 


d/n 


Sample 














Size 


5.0 


1*. 931*0 


1*. 971*2 


5.3U30 


78.801*5 


.016 


5000 




5.01*62 


5.0975 


U.9371 


80.2735 


.Oil* 






5.0108 


U.93U5 


U.6862 


71*. 9158 


.005 






U.9912 


5.0533 


5.1257 


80.9381* 


.001* 






5.0038 


5.01*12 


5.7602 


87.0035 


.003 






5.0126 


5.0367 


5.339U 


83.8973 


.005 






5.031*6 


1*. 981*1* 


U.8206 


80.1721* 


.010 






5.0016 


5.01U2 


5.1778 


77.761*8 


.006 






5.0211* 


U.9U19 


U.7638 


75.7139 


.006 






H.9598 


5.oiliO 


U.8610 


77.9821* 


.009 




Theor. 


5.0 


5.o 


5.0 


80.0 


.017 


(1050 




9.991*8 


10.0720 


10.U326 


303.7092 


.006 


5000 




10.01*02 


9.9338 


9.U176 


292 .8100 


.011 






10.0171* 


10.1211 


12.321*6 


335.7205 


.001* 






10.0260 


9.8721 


9.3898 


296.6628 


.016 






10.0852 


9.85U7 


9.8097 


302.7969 


.017 






10.0176 


9.8136 


9.8867 


29U.7293 


.011 






9.991U 


10.1918 


11.8066 


323. 7261* 


.006 




Theor. 


10.0 


10.0 


10.0 


310.0 


.017 


(1056) 
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Appendix XI. 

Exponential random number generator. 

The numbers are generated by a call to 'EXPRN(DUM) 1 . The argument 
'DUM' is not used. For example: Y * EXPRN(DUM) + U.2 
No external calls are made. The observed average generation time 
for one number is 2270 microseconds. A conversion for numbers 
with parameter other than one is given on page 31 in Section 8.1. 
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Appendix XII. 

Test results for the exponential generator. 

The test consists of consecutive samples of 1000 each for lambda 



equal one. 



SAMPLE 


Ml 


M2 


M3 


MU 


d/n Generation 

Time 


1000-1 


1.0056 


.95 UU 


1.5731 


6.0113 


.015 


-2 


1.0602 


1.1981 


2.93UO 


15.0U79 


.020 


-3 


1.0606 


I.2050 


2.9526 


15.5327 


.037 


-U 


1.0272 


1.1211 


2.5961 


12 .1087 


.020 




.9U22 


.8821 


1.7657 


7.6101 


.030 


-6 


1.02UU 


.9U8U 


1.U382 


5.209U 


.031 


-7 


.9273 


.8U57 


1.U637 


5.7615 


.038 


-8 


1.0520 


1.1072 


2.U988 


12.2010 


.033 


-9 


1.0323 


1.2839 


3.8978 


23.652U 


.016 


-10 


.9751 


.8368 


1.3691 


5.3359 


.022 2270 microsecs 


Theor. 


1.0 


1.0 


2.0 


9.0 


.039 (1056) 


-11 


.985U 


1.0562 


2.U39U 


11.612U 


.018 


-12 


.9998 


1.0U95 


2.2U30 


10.3795 


.01U 


-13 


.98U1 


1.0UU7 


2.2683 


10.1315 


.037 


-1U 


1.020U 


1.0517 


1.9077 


7.1153 


.021 


-15 


.9513 


.7981 


1.186U 


U.1650 


.022 


-16 


.9592 


.9722 


1.7937 


6.9373 


.0U1 


-17 


1.0U03 


1.1069 


2.2U5U 


9.935U 


.017 


-18 


1.0329 


1.0635 


1.9U97 


7.6527 


.018 


-19 


.9711 


.901U 


1.6697 


6.832U 


.022 


-20 


.9835 


.9U51 


1.7UU5 


7.1993 


.020 


-21 


1.0051 


1.0281 


2.15U7 


9.U030 


.015 


-22 


.9631 


.8759 


1.3956 


U.915U 


.020 


-23 


.9951 


.9813 


1.778U 


7.1579 


.015 


-2U 


1.0176 


1.1237 


2.5168 


12 .1998 


.020 


-25 


1.0278 


1.066U 


2.2520 


10.2068 


.025 


-26 


1.0317 


1.0396 


2.1303 


9.3136 


.027 


-27 


1.0188 


.9595 


1.6659 


6.7573 


.02U 


-28 


.933U 


.8338 


1.3261 


U.7659 


.038 


-29 


.9 U68 


.9008 


1.8U11 


8.U889 


.03U 


-30 


1.0172 


.9759 


1.8om 


7.57U6 


.029 
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Addendum 1 



Results of the investigation of the tails of the Marsaglia normal 
random number generator. 

The test results belie the idea resulting from investigations 
made in Section U.3* A typical series of results of the tests are 
shown below. The first line gives the theoretical cumulative sample 
result based on a sample size of 100, The following data shows 
the experimental results. The first interval is from minus infinity 
to -2.56, following intervals are 0.16 in width. Thus only the 



negative half of the distribution is Shown. 
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